﻿using System;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Decomposition;
using NUnit.Framework;

namespace dnAnalytics.Tests.LinearAlgebra.Decomposition
{
    public abstract class SvdTests
    {
        private const string singularMatrixFile = "./TestData/Matrices/gear_integer_general_coordinate_100.mtx";
        private const string squareMatrixFile = "./TestData/Matrices/random_real_general_array_100.mtx";
        private const string tallMatrixFile = "./TestData/Matrices/random_real_general_array_20_10.mtx";
        private const string wideMatrixFile = "./TestData/Matrices/random_real_general_array_10_20.mtx";

        private Matrix singularMatrix;
        private Svd singularSvd;

        private Matrix squareMatrix;
        private Svd squareSvd;
        private Matrix tallMatrix;
        private Svd tallSvd;
        private Matrix wideMatrix;
        private Svd wideSvd;

        public abstract Vector GetVector(int size);
        public abstract Matrix GetMatrix(int order);
        public abstract Matrix GetMatrix(int rows, int columns);
        public abstract Matrix GetMatrix(string matrixFile);

        [TestFixtureSetUp]
        public void SetUp()
        {
            squareMatrix = GetMatrix(squareMatrixFile);
            singularMatrix = GetMatrix(singularMatrixFile);
            wideMatrix = GetMatrix(wideMatrixFile);
            tallMatrix = GetMatrix(tallMatrixFile);

            squareSvd = new Svd(squareMatrix, true);
            singularSvd = new Svd(singularMatrix, true);
            wideSvd = new Svd(wideMatrix, true);
            tallSvd = new Svd(tallMatrix, true);
        }

        [Test]
        [ExpectedException(typeof (ArgumentNullException))]
        public void SvdConstructorNull()
        {
            new Svd(null, false);
        }

        [Test]
        public void Decomposition()
        {
            Matrix u = squareSvd.U();
            Matrix w = squareSvd.W();
            Matrix vt = squareSvd.VT();
            Matrix result = u*w*vt;

            double error;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = System.Math.Abs(result[i, j] - squareMatrix[i, j]);
                    Assert.IsTrue(error < Constants.AcceptableError);
                }
            }

            u = tallSvd.U();
            w = tallSvd.W();
            vt = tallSvd.VT();
            result = u*w*vt;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = System.Math.Abs(result[i, j] - tallMatrix[i, j]);
                    Assert.IsTrue(error < Constants.AcceptableError);
                }
            }

            u = wideSvd.U();
            w = wideSvd.W();
            vt = wideSvd.VT();
            result = u*w*vt;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = System.Math.Abs(result[i, j] - wideMatrix[i, j]);
                    Assert.IsTrue(error < Constants.AcceptableError);
                }
            }

            u = singularSvd.U();
            w = singularSvd.W();
            vt = singularSvd.VT();
            result = u*w*vt;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = System.Math.Abs(result[i, j] - singularMatrix[i, j]);
                    Assert.IsTrue(error < Constants.AcceptableError);
                }
            }
        }


        [Test]
        public void Rank()
        {
            int rank = squareSvd.Rank();
            Assert.AreEqual(squareMatrix.Rows, rank);

            rank = tallSvd.Rank();
            Assert.AreEqual(tallMatrix.Columns, rank);

            rank = wideSvd.Rank();
            Assert.AreEqual(wideMatrix.Rows, rank);

            rank = singularSvd.Rank();
            Assert.AreEqual(squareMatrix.Rows - 1, rank);
        }

        [Test]
        public void Norm()
        {
            double norm = squareSvd.Norm2();
            Assert.AreEqual(1011.4190766097751, norm, Constants.AcceptableError);

            norm = tallSvd.Norm2();
            Assert.AreEqual(445.27398550522514, norm, Constants.AcceptableError);

            norm = wideSvd.Norm2();
            Assert.AreEqual(385.51714390506197, norm, Constants.AcceptableError);

            norm = singularSvd.Norm2();
            Assert.AreEqual(2, norm, Constants.AcceptableError);
        }

        [Test]
        public void ConditionNumber()
        {
            double cn = squareSvd.ConditionNumber();
            Assert.AreEqual(207.5103850393096, cn, 1e-4);

            cn = tallSvd.ConditionNumber();
            Assert.AreEqual(5.01902319723716, cn, Constants.AcceptableError);

            cn = wideSvd.ConditionNumber();
            Assert.AreEqual(3.65923703921971, cn, Constants.AcceptableError);

            cn = singularSvd.ConditionNumber();
            Assert.Greater(cn, 1e16);
        }
    }
}